function estimate_pf_TL(treatment,k) 

%%% Settings
nMinObs = 25; % minimum number of obs required for GMM
dateBeginSample = 140; %corresponds to 1995Q1

% backward-looking rolling window for data (no. of quarters incl. current)
movingWindowSize = 20; 

% optimizer options
optionsUnc = optimoptions(@fminunc,'Display','off','StepTolerance',1e-14,'FunctionTolerance',1e-14); 
optionsNM = optimset('Display','off'); 

%%% load data
dataIn = readtable(['../../B_temp/pf_input_naicsq' num2str(k) '_' treatment '.csv']);

% create lists of industries and dates  
allNaicsCodes   = unique(dataIn{:,'naicsq'});   
allDates        = unique(dataIn{:,'date'}); 

%%% estimation

N = length(allNaicsCodes); 
for i = 1:N 

        % industry-quarter considered
        thisNaicsCode = allNaicsCodes(i); 

        % save estimates: indices
        naicsList(i,1) = thisNaicsCode; 

        % select data
        relevantRows = (dataIn.naicsq == thisNaicsCode & dataIn.date >= dateBeginSample ); 

        % start estimation if observations available   
        if sum(relevantRows) < nMinObs

            betaCogsList(i,1) = NaN;
            betaCapitalList(i,1) = NaN;
            betaCogs2List(i,1) = NaN;
            betaCapital2List(i,1) = NaN;
            thetaV(i,1)=NaN;
            thetaK(i,1)=NaN;
            
        else
            
            % data and OLS estimate
            data = dataIn{relevantRows,{'y','v','k'}};    
                
            % first-stage measurement error correction
            dataPoly = ones(size(data(:,2)));
            for ii=1:3
                dataPoly = [dataPoly data(:,2).^ii];
                dataPoly = [dataPoly data(:,3).^ii];
                for jj=1:3
                    dataPoly = [dataPoly data(:,2).^ii .* data(:,3).^jj];
                end
            end  
                
            coeff = dataPoly\data(:,1); 
            y_corr = dataPoly*coeff;
            Ly_corr = [NaN; y_corr(1:end-1)];

            % create data with lags in panel structure
            id = dataIn{relevantRows,{'id'}}; 
            tt = dataIn{relevantRows,{'date'}}; 

            data = [y_corr, ...
                    dataIn{relevantRows,{'v','k'}}, ...
                    Ly_corr,...
                    dataIn{relevantRows,{'Lv','Lk'}},...  
                    dataIn{relevantRows,{'L4v'}},... 
                    id,...
                    [NaN; id(1:end-1)],...
                    tt,...
                    [NaN; tt(1:end-1)],...
                    ];

            data = data(data(:,end-1)-1 == data(:,end) & data(:,end-3) == data(:,end-2),1:end-4); 
            
            % OLS
            paramsStartOLS = ols(data(:,1),[data(:,2), data(:,3), data(:,2).^2, data(:,3).^2]); 
            paramsStartOLS = paramsStartOLS(1:end-1);

            % GMM
            [paramsUnc, fvalUnc, flagUnc] = fminunc(@(params) GMM_objective_TL(params, data), paramsStartOLS, optionsUnc);            
            [paramsNM, fvalNM, flagNM]  = fminsearch(@(params) GMM_objective_TL(params, data), paramsStartOLS, optionsNM);

            % save estimates  
            if fvalUnc <= min([fvalUnc,fvalNM])
                betaCogsList(i,1) = paramsUnc(1);
                betaCapitalList(i,1) = paramsUnc(2); 
                betaCogs2List(i,1) = paramsUnc(3); 
                betaCapital2List(i,1) = paramsUnc(4); 
                method = 'unc';
            elseif fvalNM <= min([fvalUnc,fvalNM])
                betaCogsList(i,1) = paramsNM(1);
                betaCapitalList(i,1) = paramsNM(2); 
                betaCogs2List(i,1) = paramsNM(3); 
                betaCapital2List(i,1) = paramsNM(4); 
                method = 'NM';
            end 

    end 
    
end

%% results
results = table(naicsList,betaCogsList,betaCogs2List);
results.Properties.VariableNames = {'naics', 'betaCogsTL','betaCogs2TL'}; 
writetable(results,['../../B_temp/pf_TL_naicsq' num2str(k) '_' treatment '.csv'],'Delimiter',',')

end